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Abstract 

We analyze the stability of synchronized state for coupled nearly identical dynamical systems on 
networks by deriving an approximate Master Stability Function (MSF). Using this MSF we treat 
the problem of designing a network having the best synchronizability properties. We find that the 
edges which connect nodes with a larger relative parameter mismatch are preferred and the nodes 
having values at one extreme of the parameter mismatch are preferred as hubs. 
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When two or more dynamical systems are coupled or driven by a common signal the 



systems may synchronize under suitable conditions 
synchronization, such as complete synchronization 



ll, |2|. One can achieve different types of 

n n 

2j, phase synchronization [3[, lag syn- 



chronization [4|, generalized synchronization js] etc. Recently, there is considerable interest 
in the synchronization of coupled dynamical systems on a network [6| . For coupled identical 
systems which give exact synchronization, Pecora and Carroll j3] have introduced a mas- 
ter stability function (MSF) which can be calculated from a simple set of master stability 
equations and then applied to the study of stability of the synchronous state of different 
networks. This general approach has become popular and has been used in various studies 



of synchronization on networks 8l4l2l| . Several works on different networks have shown that 
small world and scale free networks show better synchronization properties js, 9|. 

For coupled nonidentical systems, in general it is difficult to obtain exact synchronization. 
But, one can get synchronization of some generalized type [5]. The parameter mismatch 
between different coupled systems can lead to desynchronization bursts and this is known 
as the bubbling transition [l3, [l^. Restrepo et el. 15| have studied the spatial patterns of 



such desynchronization bursts in networks. After the desynchronization burst the system 



returns to the synchronized state. Sun et al. [16| determine the deviation from average 
trajectory as a function of the mismatch. However, for nonidentical systems, there is no 
general theory such as MSF, to study the stability of synchronization. 

In this paper we address the question of the stability of synchronization of coupled nearly 
identical systems on networks. By using the property of differential equations that the ho- 
mogeneous part determines the exponential rates and treating the parameter mismatch in a 
first order perturbation theory, we derive master stability equation for coupled nearly identi- 
cal systems. This master stability equation uses the homogeneous state and two parameters, 
a for the network coupling and A for the mismatch in nonidentical systems. This allows us 
to define the MSF and study the stability properties of the synchronized state. 

When one considers identical coupled systems the important question is about the type of 
network which gives better synchronization properties. When one considers coupled nearly 
identical systems, additional interesting and important questions arise. Which nodes are 
better chosen as hubs? Which edges give better synchronization? Using our MSF we find 
that for better synchronization nodes on one extreme of parameter mismatch are preferred 
as hubs and nodes with larger relative parameter mismatch are preferred for constructing 



edges. 

Consider coupled dynamical systems, 

N 

x' = f{x\ rO + £ 5^ G,,h{x^); t = l,...,N (1) 
j=i 

where, x*(g R"^) is an m dimensional state vector of the system i, f : R"^ — )■ R^ gives 
the dynamics of an isolated system, e is a scalar coupling parameter and h : R"^ — )■ i?"^ is 
a coupling function, G is the coupling matrix of the network, is some parameter which 
depends on the node i. 

For the coupled identical systems, i.e. = r, V?, the synchronization manifold is defined 
by = ■ ■ ■ = x^ = X and is an invariant manifold provided the coupling matrix satisfies the 
condition that Gij = 0, Vi. With this condition, the synchronized state x, is a solution 
of the uncoupled dynamics, x = f{x). 

The condition J^jGij = ensures that G has one eigenvector ei = (1,...,!)-^, with 
eigenvalue 71 = 0. This eigenvector defines the synchronization manifold. All the remaining 
eigenvectors belong to the transverse manifold. The synchronized state is stable provided 
all the transverse Lyapunov exponents are negative. 

Now, let us consider the case when the parameter depends on the node i. Let the 
parameter mismatch be Sr^ = ri — r where r is some typical value of the parameters r,. 
In general, for nonidentical systems it is not possible to get an exact synchronization of 
the type discussed above. Instead we get a generalized synchronization where there is a 
functional relationship between variables of the systems, e.g. g{x^,x^) = 0. The generalized 
synchronization is stable provided the largest transverse Lyapunov exponent is negative. 

To determining the stability of this generalized synchronization, we do the linear stability 
analysis. In this analysis, we retain terms to second order in 2;* = — x and 5rj. The reason 
for doing this will be clear shorty. Thus the dynamics of the deviation can be written as 

= D^f{x, f)z' + GijD^h{x)z^ + Drf{x, f)5ri 

+DrDJ{x, r)z'6n + ]^Dlf{x, f)5r} + . . . (2) 

The terms corresponding to {z^Y included since we will be interested in the solution 

2;* = for finite (5rj. As an equation for 2;*, the RHS of Eq. (|2]) contains both homogeneous 
and inhomogeneous terms. To a first approximation, the inhomogeneity won't affect the 
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exponential rate of convergence of the trajectories to the synchronous solutions though it 
can shift the solution. To see this consider a general linear equation Du = p{t), where D 
is a differential operator. Let the solution he u = Uh + g(t) where Uh = Aihi{t)exp{kit) 
is the solution of the homogeneous equation Du = 0, and Ai are constants. If p{t) does 
not have any exponential dependence, then g(t) cannot contain any additional exponential 
other than already in Uh, due to the property that the derivative of an exponential is also 
an exponential with the same exponent. For example, for ii = —ku + p, the solution of the 
homogeneous equation is Uh{t) = u{0)e~^^ and of the inhomogeneous equation with constant 
p is u{t) = (u(0) — {p/k))e^^* +p/k. We note that the inhomogeneity in the differential 
equation shifts the asymptotic solution but does not change the exponential. In our case the 
stability of the synchronized state is governed by the largest transverse Lyapunov exponent, 
i.e. only by the exponential rates which are determined by the homogeneous equation. The 
inhomogeneous part will shift the solution. In addition, while calculating the Lyapunov 
exponents, it is necessary that the shifted solution preserves the nature of the attractor so 
that the average expansion and contraction rates are not significantly affected by the shift. 
This can be assumed to be valid when different systems are in generalized synchrony since 
they are related to each other. This may also hold very near the synchronization region but 
not far away from it. 

Hence, to obtain the Lyapunov exponents, we consider the homogeneous equation ob- 
tained from Eq. (jj]). 

N 

= DJ z' + eY, G.,jD,h z^ + DrDJ z'5n (3) 

This equation can be put in a matrix form as 

Z = DJ Z + sD^h Z + DrD^f Z R (4) 

where Z = {z^, . . . , z^) is an m x iV matrix and R = diag((5ri, . . . , 6t^) is an N x N diagonal 
matrix. 

Let 7fc, , = 1, . . . , be the eigenvalues and right eigenvectors of G^. Acting Eq. (jl]) 
on and using the m dimensional vectors (pk = Ze§, we get 

k = [DJ + eikDM<Pk + DrDJ Z R ef . (5) 
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In general, are not eigenvectors of R and hence Eq. ([5]) is not easy to treat. To solve 
Eq. ([5]) we use first order perturbation theory and write Eq. ([5]) as 

0fc = [Dxf + e-fkD^u + UkDrD^f](j)k (6) 

where = e^Re^ is the first order correction and is the left eigenvector of G^. 

Since both 7^ and z/^ can be complex, treating them as complex parameters a = e'jk and 
A = z/fc respectively, we can construct the master stability equation as 

= [DJ + aD^h + ADrDJ]<j). (7) 

For the coupled identical systems, the above equation reduces to the master stability equa- 
tion given by Pecora and Carroll [7]. We can determine the MSF or Amax, which is the 
lar gest Lyapunov exponent for Eq. ((Tj), as a surface in the complex space defined by a and 



A [2l|]. The synchronized state is stable if the MSF is negative at each of the eigenvalues 
7fc = a and = A {k ^ 1). This ensures that all the transverse Lyapunov exponents are 
negative. 

We note that though the master stability equation uses the homogeneous state, it 
allows us to study the stability of the generalized synchronization in nonidentical systems. 
The mismatch between the different systems is included through the parameter A. 

To examine how well Eq. ([7]) allows the estimation of Lyapunov exponents, we calcu- 
late the Lyapunov exponents for the coupled Rossler systems |20| and compare them with 
those obtained from Eq. ([7]). Consider coupled chaotic Rossler systems with different 
frequencies, 

N 

iji = ujiXi + a^Vi (8) 

where Ui is the frequencies of the i-th oscillator and Lij = 1 if the nodes i and j are coupled 
and zero otherwise and La = — ^j^iLij. For simplicity we restrict ourselves to symmetric 
coupling matrices L so that the eigenvalues and hence a and A are real. 

We first consider two coupled Rossler oscillators. Fig. 1 plots the three largest Lyapunov 
exponents, Aj, 2 = 1,2, 3, as a function of the coupling strength e and their estimated values 
fj^Qjji Eq. ([7]). Fig. 2a plots the difference 5\i = Aj — A^^"^ as a function of e for these 
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FIG. 1. The figure shows the three largest Lyapunov exponents Aj, z = 1, 2, 3 (red, green and blue) 
and their estimated values obtained from the master stability equation (Eq. ([7])) (pink, cyan 
and black) as a function of e for two coupled Rossler systems with frequencies uji = 1.05 and 
UJ2 = 1.07. Taking Co = 1.0 we get Ai = A2 = 0.06 which are used in Eq. d?]). Rossler parameters 
are Ur = bj. = 0.2, Cr = 7.0. The synchronous state is stable in the region given by ai < a < 02 
indicated by the arrows. 

Lyapunov exponents. The region when the third largest Lyapunov exponent A3 < 0, corre- 
sponds to the synchronization region and in this region it is the largest transverse Lyapunov 
exponent. From Figs. 2a, we find that the differences 6Xi are small in the synchronization 
region and very close to it. Though only three exponents are plotted in the figure, the 
differences are small for the other Lyapunov exponents. Fig. 2b plots the difference 6Xi as 
a function of e for the three largest Lyapunov exponents for a random network of sixteen 
nodes. Again we observe that the errors are small in the synchronization region. Thus, we 
find that the master stability equation ([7]) can estimate the actual Lyapunov exponents for 
the synchronized state reasonably well. 

Now, we consider the MSF, Xmax, which is the largest transverse Lyapunov exponent. 
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FIG. 2. a. The figure shows the difference 5\i = Xi — for the three largest Lyapunov exponents 
as a function of the couphng constant e for two coupled Rossler systems with parameters as in 
Fig 1. b. The figure shows the difference 6Xi for the three largest Lyapunov exponents as a function 
of e for sixteen randomly coupled Rossler systems having different internal frequencies uji . We find 
that the differences are small in the synchronization region. 

It can be calculated using Eq. ([7]). In Fig. [3] we plot Xmax in the parameter plane (a, A) 
as a contour plot for Rossler system. From the figure we can see that the stability region 
increases with the parameter A. 

We now demonstrate the utility of the master stability function by considering the prob- 
lem of construction of an optimized network which gives best synchronization properties. 
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FIG. 3. The master stability function Xmax for Rossler system is plotted as a contour plot in the 
parameter plane (a, A). The stability region is given by the "V" shape region bordered by the 
contours from both sides. 



To construct the optimized network we adapt Monte Carlo optimization method [22|] and 
rewire the edges of the network to construct a network that shows best synchronizability, 
i.e. the largest interval 1^ of the coupling constant e which shows synchronization. 

We start with a system of nearly identical coupled Rossler oscillators as in Eq. on 
a connected network of nodes and E randomly chosen edges. In each Monte Carlo 
step we rewire one edge. If the rewired network increases the stability interval 1^ of the 
synchronized state, then it is chosen with probability one, otherwise it is accepted with 

8 



probability ^PiiT"" -^f"^) where (3 is the inverse temperature. 

We now investigate two questions. In the optimized network, which edges are more 
preferable and which nodes have larger number of connection or act as hubs? 

To investigate the question of which nodes act as hubs, we define the correlation coefficient 
between the frequency and the degree of a node as pi^k = <(fc,-<fc,>)(cj,-<^,>)> where 
ki = —La is the degree of node i. Fig. shows p^jfe (solid line) as a function of Monte 
Carlo steps. For the random network p^^ = 0. We find that p^^ increases and saturates to 
a positive value. Thus, in the synchronized optimized network the nodes which have larger 
frequencies have more connections and are preferred as hubs. The reason for this is the 
"V" shape of the stability region in Fig. [3l i.e. the stability range increases as A increases. 
We have also investigated a case were an opposite behavior is obtained. If instead of the 
frequency, we make the parameter in Eq. (|8]) node dependent, then the stability region in 
the plot of MSF similar to Fig. |3l has an inverted "V" shape. In this case in the optimized 
network, nodes which have smaller values of have more connections and are preferred as 
hubs. 

To investigate the question of which edges are preferred, we define the correlation co- 
efficient between the absolute frequency differences between two nodes and the edges as 
Poja = <(^'j <A,j>)(\u), u)j\ <\u), u)j\>)> where Aij = 1 if nodes i and ?' are connected and 
otherwise. Fig. HJa shows p^ia as a function of Monte Carlo steps. We find that p^a in- 
creases from (the value for the random network) and saturates. Thus, in the synchronized 
optimized network the pair of nodes which have a larger relative frequency mismatch are 
preferred as edges for the optimized network. Again, the reason for this preference of edges 
is probably the conical shape of the stability region in Fig. 3. The edges are to be chosen 
so that the parameter A increases and the stability region increases. 

To conclude we have developed the Master Stability Function (MSF) approach for coupled 
nonidentical systems. We use the property of differential equations that the homogeneous 
part is mainly responsible for the exponential dependence of the variables. The parameter 
mismatch is treated in a first order perturbation theory. Our MSF uses the homogeneous 
state but it still allows us to study the stability properties of generalized synchronization 
for nonidentical systems. Using MSF, we construct optimized networks with better syn- 
chronization properties by rewiring the network keeping the number of edges constant. We 
find that in the optimized network the nodes having parameter mismatch at one extreme 
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FIG. 4. The figure plots the correlation coefficient p^^k (figure a) and pu)a (figure b) as a function 
of the Monte Carlo steps of optimization for 32 coupled Rossler systems. We see that both p^^ 
and Pijja increase and saturate to positive values. 

depending on the shape of stabihty region in MSF plot, have more edges and are preferred 
as hubs and the pair of nodes which have a larger relative parameter mismatch are preferred 
for constructing edges. 
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